- Plot vector geospatial data (points/polygons)
- Plot raster geospatial data (image)
- Extracting data from raster
- Making maps (Compass rose, overlays, etc)
2020-02-19
We are researchers who want to answer the question "What range of environmental conditions do taxon tend to inhabit?"
General workflow:
Extract data about environmental conditions at those locations
Variants of this question/workflow are applicable to others study systems.
–> –> –> –> –> –> –>
–> –> –> –> –> –> –> –> –> –> –> –> –>
rgdal released in 2003–import from more geographic data formatssp released in 2005–creates spatial objects with classes and generic methods for points, lines, polygons, grids, and attribute data (but still lacked ability to do geometric operations)raster released in 2010–support for raster operations and moreGRASS, spgrass6, and rgrass7RSAGA, RPyGeo, and RQGISggplot2ggmaprasterVistmapleafletmapview sfsp, rdgal for read/write, and rgeos for geometric operationstidyversesf packagesf (simple feature) objects are extended data.frames or tibblestibblesfc (simple feature columns) with bounding box bbox and CRS (coordinate reference system) attributes## Reading layer `NAM-level_1' from data source `/Users/michelebuonanduci/Desktop/QERM597-making-maps/data/NAM-level_1.shp' using driver `ESRI Shapefile' ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : ## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed ## incompletely to real 0. ## Simple feature collection with 19 features and 4 fields ## geometry type: POLYGON ## dimension: XY ## bbox: xmin: 11.71563 ymin: -28.96946 xmax: 25.2567 ymax: -16.95989 ## epsg (SRID): NA ## proj4string: NA ## [1] "sf" "data.frame" ## Simple feature collection with 6 features and 4 fields ## geometry type: POLYGON ## dimension: XY ## bbox: xmin: 13.62839 ymin: -28.96946 xmax: 25.2567 ymax: -17.47639 ## epsg (SRID): NA ## proj4string: NA ## AREA PERIMETER ID CAPTION geometry ## 1 NA 0 Caprivi Caprivi POLYGON ((24.77467 -17.8854... ## 2 NA 0 Erongo Erongo POLYGON ((15.81059 -23.2945... ## 3 NA 0 Hardap Hardap POLYGON ((18.48226 -25.4376... ## 4 NA 0 Karas Karas POLYGON ((16.3988 -25.64753... ## 5 NA 0 Karas <NA> POLYGON ((15.20197 -27.0306... ## 6 NA 0 Karas <NA> POLYGON ((14.92512 -26.3011... ## Reading layer `enproads' from data source `/Users/michelebuonanduci/Desktop/QERM597-making-maps/data/enproads.shp' using driver `ESRI Shapefile' ## Simple feature collection with 565 features and 3 fields ## geometry type: LINESTRING ## dimension: XY ## bbox: xmin: 14.35337 ymin: -19.48418 xmax: 17.15124 ymax: -18.5008 ## epsg (SRID): NA ## proj4string: NA ## [1] "sf" "data.frame" ## Simple feature collection with 6 features and 3 fields ## geometry type: LINESTRING ## dimension: XY ## bbox: xmin: 14.46841 ymin: -19.41752 xmax: 14.58964 ymax: -19.33077 ## epsg (SRID): NA ## proj4string: NA ## TYPE LENGTH KM geometry ## 1 Track 0.7 1 LINESTRING (14.50502 -19.35... ## 2 Track 3.6 4 LINESTRING (14.49077 -19.36... ## 3 Graded 2.3 2 LINESTRING (14.4864 -19.417... ## 4 Track 7.1 7 LINESTRING (14.58498 -19.34... ## 5 Track 1.2 1 LINESTRING (14.58498 -19.34... ## 6 Graded 0.1 0 LINESTRING (14.52596 -19.33...
sf and tidyversesf functions begin with st_summary, plotsf methods for filter, arrange, distinct, group_by, ungroup, mutatemethods(class = "sf") ## [1] [ [[<- ## [3] $<- aggregate ## [5] anti_join arrange ## [7] as.data.frame cbind ## [9] coerce dbDataType ## [11] dbWriteTable df_spatial ## [13] distinct extent ## [15] extract filter ## [17] full_join gather ## [19] group_by group_map ## [21] group_split identify ## [23] initialize inner_join ## [25] left_join mask ## [27] merge mutate ## [29] nest plot ## [31] print raster ## [33] rasterize rbind ## [35] rename right_join ## [37] sample_frac sample_n ## [39] select semi_join ## [41] separate show ## [43] slice slotsFromS3 ## [45] spread st_agr ## [47] st_agr<- st_area ## [49] st_as_sf st_bbox ## [51] st_boundary st_buffer ## [53] st_cast st_centroid ## [55] st_collection_extract st_convex_hull ## [57] st_coordinates st_crop ## [59] st_crs st_crs<- ## [61] st_difference st_force_polygon_cw ## [63] st_geometry st_geometry<- ## [65] st_interpolate_aw st_intersection ## [67] st_intersects st_is_polygon_cw ## [69] st_is st_line_merge ## [71] st_linesubstring st_make_valid ## [73] st_minimum_bounding_circle st_nearest_points ## [75] st_node st_normalize ## [77] st_point_on_surface st_polygonize ## [79] st_precision st_segmentize ## [81] st_set_precision st_simplify ## [83] st_snap_to_grid st_snap ## [85] st_split st_subdivide ## [87] st_sym_difference st_transform_proj ## [89] st_transform st_triangulate ## [91] st_union st_voronoi ## [93] st_wrap_dateline st_write ## [95] st_zm summarise ## [97] transmute ungroup ## [99] unite unnest ## see '?methods' for accessing help and source code
st_geometry(Namibia) %>% plot
plot(st_geometry(roads))
# Default methods for objects plot(Namibia) ## Warning in min(x): no non-missing arguments to min; returning Inf ## Warning in max(x): no non-missing arguments to max; returning -Inf
plot(Namibia[,3])
plot(roads)
plot(roads[,1])
plot(filter(roads[,3], KM > 20))
fileloc <- system.file("shape/nc.shp", package = "sf")
nc <- read_sf(fileloc)
st_crs("+proj=longlat +datum=WGS84") # Proj.4 string"
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(3857) # ESPG code: public registry of spatial reference systems, 3857 = web-based mapping (i.e. Google maps, OpenStreetMap)
## Coordinate Reference System:
## EPSG: 3857
## proj4string: "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
st_crs(3857)$units # check units
## [1] "m"
st_crs(NA) # unknown (assumed planar/Cartesian)
## Coordinate Reference System: NA
units::set_units(a1, km^2) ## 1137.389 [km^2] units::set_units(a2, km^2) ## 1137.577 [km^2] units::set_units(a3, km^2) ## 1137.577 [km^2]
Spatial* objects using st_as_sfrnaturalearthhead(hotsprings_df) ## # A tibble: 6 x 8 ## STATE LAT LONG SpringName Temp_F Temp_C AREA USGS_quad ## <chr> <dbl> <dbl> <chr> <dbl> <dbl> <chr> <chr> ## 1 CA 38.8 -123. LITTLE GEYS… 210 99 SANTA … (WHISPERING … ## 2 CA 41.5 -120. HOT SPRINGS… 208 98 ALTURAS CEDARVILLE 15 ## 3 CA 36.0 -118. COSO HOT SP… 207 97 DEATH … HAIWEE RESER… ## 4 CA 41.7 -120. LAKE CITY H… 207 97 ALTURAS CEDARVILLE 15 ## 5 CA 36.0 -118. DEVILS KITC… 207 97 DEATH … HAIWEE RESER… ## 6 CA 40.4 -121. TERMINAL GE… 205 96 SUSANV… MT. HARKNESS…
crs = 4326)hotsprings <- st_as_sf(hotsprings_df,
coords = c("LONG", "LAT"),
crs = 4326)
head(hotsprings)
## Simple feature collection with 6 features and 6 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -122.748 ymin: 36.036 xmax: -117.769 ymax: 41.67
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 7
## STATE SpringName Temp_F Temp_C AREA USGS_quad geometry
## <chr> <chr> <dbl> <dbl> <chr> <chr> <POINT [°]>
## 1 CA LITTLE … 210 99 SANTA… (WHISPER… (-122.748 38.767)
## 2 CA HOT … 208 98 ALTUR… CEDARVIL… (-120.078 41.534)
## 3 CA COSO … 207 97 DEATH… HAIWEE R… (-117.769 36.045)
## 4 CA LAKE … 207 97 ALTUR… CEDARVIL… (-120.206 41.67)
## 5 CA DEVILS … 207 97 DEATH… HAIWEE R… (-117.802 36.036)
## 6 CA TERMINAL … 205 96 SUSAN… MT. HARK… (-121.375 40.421)
world <- ne_countries(scale = "medium", returnclass = "sf") st_geometry(world) %>% plot()
st_writeZebra.csv and plot the geometry.roads <- st_read(here("data", "enproads.shp"), crs = "+init=epsg:4326") %>% #4326= coordinate system based on Earth's center of mass
st_transform("+init=epsg:32733") #32733 = spatial reference for Nambia
## Reading layer `enproads' from data source `/Users/michelebuonanduci/Desktop/QERM597-making-maps/data/enproads.shp' using driver `ESRI Shapefile'
## Simple feature collection with 565 features and 3 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 14.35337 ymin: -19.48418 xmax: 17.15124 ymax: -18.5008
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
st_geometry(roads) %>% plot
zebra_sf <- read.csv(here("data", "Zebra.csv")) %>%
dplyr::select(ID = Name, 4:6) %>%
mutate(timestamp = as.POSIXct(lubridate::mdy_hm(Date))) %>%
st_as_sf(., coords = 3:4, crs = "+init=epsg:4326") %>% #longlat, converting foreign object to SF object
st_transform("+init=epsg:32733") #convert to UTM
ggplot() +
geom_sf(data=roads) +
geom_sf(data=zebra_sf, aes(color=ID))
How close or far from different road types do zebra move?
unique(roads$TYPE)
## [1] Track Graded Gravel Tar
## Levels: Graded Gravel Tar Track
#Filter roads based off type:
large_roads <- filter(roads, TYPE %in% c("Tar", "Gravel"))
small_roads <- filter(roads, TYPE %in% c("Graded", "Track"))
ggplot() +
geom_sf(data=large_roads, size=1.5) +
geom_sf(data=small_roads, size=0.6) +
geom_sf(data=zebra_sf, aes(color=ID))
# Find the minimum distance (in meters) of each point to a large road
large_dist<- st_distance(y=zebra_sf, x=large_roads) # a units matrix dim =[nrow(x), nrow(y)]; takes 10-20 seconds
zebra_sf$large_road_dist <- apply(large_dist, 2, min)
# Find the minimum distance (in meters) of each point to a small road
small_dist<- st_distance(y=zebra_sf, x=small_roads) # a units matrix dim =[nrow(x), nrow(y)]; takes 10-20 seconds
zebra_sf$small_road_dist <- apply(small_dist, 2, min)
head(data.frame(zebra_sf))
## ID Date timestamp geometry
## 1 AG059 4/20/2009 15:40 2009-04-20 15:40:00 POINT (594243.6 7885501)
## 2 AG059 4/20/2009 13:01 2009-04-20 13:01:00 POINT (596576 7879266)
## 3 AG059 4/20/2009 12:02 2009-04-20 12:02:00 POINT (596579.5 7879260)
## 4 AG059 4/20/2009 17:00 2009-04-20 17:00:00 POINT (584733 7890124)
## 5 AG059 4/21/2009 1:40 2009-04-21 01:40:00 POINT (596611.3 7879269)
## 6 AG059 4/21/2009 1:00 2009-04-21 01:00:00 POINT (596605.8 7879266)
## large_road_dist small_road_dist
## 1 7.845860 3479.8987
## 2 163.187201 241.8660
## 3 156.746071 245.8446
## 4 9.115608 1605.2154
## 5 151.769751 227.9759
## 6 151.358102 231.6984
ggplot(zebra_sf) +
geom_histogram(aes(large_road_dist, fill="large roads"), alpha=0.5) +
geom_histogram(aes(small_road_dist, fill="small roads"), alpha=0.5) +
scale_fill_manual(labels=c("large roads", "small roads"), values=c("blue", "orange"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## class : RasterLayer ## dimensions : 474, 1274, 603876 (nrow, ncol, ncell) ## resolution : 231.6564, 231.6564 (x, y) ## extent : 431873.8, 727004, 7844489, 7954294 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : /Users/michelebuonanduci/Desktop/QERM597-making-maps/data/ndvi_mean_utm.tif ## names : ndvi_mean_utm ## values : -558.117, 5303.021 (min, max) ## class : RasterLayer ## dimensions : 474, 1274, 603876 (nrow, ncol, ncell) ## resolution : 231.6564, 231.6564 (x, y) ## extent : 431873.8, 727004, 7844489, 7954294 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : /Users/michelebuonanduci/Desktop/QERM597-making-maps/data/ndvi_mean_utm.tif ## names : ndvi_mean_utm ## values : -558.117, 5303.021 (min, max) ## class : RasterLayer ## dimensions : 474, 1274, 603876 (nrow, ncol, ncell) ## resolution : 231.6564, 231.6564 (x, y) ## extent : 431873.8, 727004, 7844489, 7954294 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : /Users/michelebuonanduci/Desktop/QERM597-making-maps/data/ndvi_mean_utm.tif ## names : ndvi_mean_utm ## values : -558.117, 5303.021 (min, max)
## class : RasterLayer ## dimensions : 474, 1274, 603876 (nrow, ncol, ncell) ## resolution : 231.6564, 231.6564 (x, y) ## extent : 431873.8, 727004, 7844489, 7954294 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : memory ## names : ndvi_mean_utm ## values : -0.0558117, 0.5303021 (min, max)
## ID Date timestamp geometry ## 1 AG059 4/20/2009 15:40 2009-04-20 15:40:00 POINT (594243.6 7885501) ## 2 AG059 4/20/2009 13:01 2009-04-20 13:01:00 POINT (596576 7879266) ## 3 AG059 4/20/2009 12:02 2009-04-20 12:02:00 POINT (596579.5 7879260) ## 4 AG059 4/20/2009 17:00 2009-04-20 17:00:00 POINT (584733 7890124) ## 5 AG059 4/21/2009 1:40 2009-04-21 01:40:00 POINT (596611.3 7879269) ## 6 AG059 4/21/2009 1:00 2009-04-21 01:00:00 POINT (596605.8 7879266) ## large_road_dist small_road_dist NDVI_mean ## 1 7.845860 3479.8987 0.2409661 ## 2 163.187201 241.8660 0.2682169 ## 3 156.746071 245.8446 0.2682169 ## 4 9.115608 1605.2154 0.2380762 ## 5 151.769751 227.9759 0.2419201 ## 6 151.358102 231.6984 0.2419201
ggmapggmap GitHub page for more information about API keysggmapggmapTwo steps:
This is my personal API key, which you can use for the purposes of this class!
register_google(key = "AIzaSyBRkw_wzDKuWZDikdD86Kmp1Sa6PzuCKFc")
To add a Stamen basemap, first define the bounding box for the basemap you would like to download.
zebra_bbox <- c(14, -20, 17.5, -18)
names(zebra_bbox) <- c("left", "bottom", "right", "top")
terr_basemap <- get_stamenmap(bbox=zebra_bbox, zoom=9, maptype="terrain") ggmap(terr_basemap) + geom_sf(data=roads, inherit.aes = FALSE) + geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) + coord_sf(crs = st_crs(4326))
wat_basemap <- get_stamenmap(bbox=zebra_bbox, zoom=9, maptype="watercolor") ggmap(wat_basemap) + geom_sf(data=roads, inherit.aes = FALSE) + geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) + coord_sf(crs = st_crs(4326))
To add a Google basemap, first define the center coordinates for the basemap you would like to download.
zebra_center <- c(15.8, -19)
names(zebra_center) <- c("lon", "lat")
sat_basemap <- get_googlemap(center=zebra_center, zoom=8, size=c(640, 640), scale=2,
maptype="satellite")
ggmap(sat_basemap) +
geom_sf(data=roads, inherit.aes = FALSE) +
geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326))
Notice that get_googlemap() returns a square basemap tile, which then sets the plotting limits of your map.
To manually set different plotting limits, pass in an empty "base layer" and set x and y limits using ggplot().
ggmap(sat_basemap, maprange=FALSE, base_layer=ggplot(aes(x=1, y=1), data=NULL)) +
xlim(14.3, 17.4) + ylim(-20, -18) + xlab("lon") + ylab("lat") +
geom_sf(data=roads, inherit.aes = FALSE) +
geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326))
RgoogleMapsDoes not seems to be compatible with ggplot2 (boo)
tif file as a raster brick, then display red/green/blue valuesOnly works well if you are mapping a large spatial extent (lacks fine resolution)
## OGR data source with driver: ESRI Shapefile ## Source: "/private/var/folders/vn/q1jnls511r35nvbxkfgqzn4h0000gn/T/RtmpqXdESv", layer: "ne_10m_rivers_lake_centerlines" ## with 1455 features ## It has 34 fields ## Integer64 fields read as strings: rivernum ne_id ## Warning in rgdal::readOGR(destdir, file_name, encoding = "UTF-8", ## stringsAsFactors = FALSE, : Dropping null geometries: 1453 ## although coordinates are longitude/latitude, st_intersection assumes that they are planar ## Warning: attribute variables are assumed to be spatially constant ## throughout all geometries